NASA-CR-l<?8flOA 


Lamar University 

BEAUMONT, TEXAS 


TECHNICAL REPORT 

( NAS A— CR— 1 98 80 A ) DAMAGE TOLERANCE N95 29369 

IN FILAMENT-WOUND GRAPHITE/EPCXY 
PRESSURE VESSELS Final Report 

(Lairar Univ.) 97 p Unc1as 


G3/37 0055367 



TEXt& 


COLLEGE OF ENGINEERING 


Damage Tolerance 
in 

Filament-Wound Graphite/Epoxy Pressure Vessels 


FINAL REPORT 
NASA Grant NAG 9-698 


July 1, 1995 


Lamar University 
College of Engineering 
Beaumont, Texas 77710 


William E. Simon, Ph.D., P.E., Principal Investigator 
Vinh D. Nguyen, Ph.D., Research Associate 
Ravi K. Chenna, Research Assistant 



TABLE OF CONTENTS 


Page 

List of Figures V1 

List of Tables viii 

Chapter 

1 . Introduction 

1.1 Low-Velocity Impact Damage in Composites 1 

1.2 Experimental Studies 2 

1.3 Analytical Studies 3 

1.3.1 Wave Front Models 3 

1.3.2 Contact Models 3 

1.4 Residual Strength Prediction Models 6 

1 . 5 Damage Characterization 8 

1.6 Outline of present work 12 

1.6.1 Organization 12 

2 . Formulations 

2.1 Statement of the problem 13 

2.2 Analysis 14 

2.2.1 Modeling Contact Behavior 15 

2.2.2 Finite Element Formulation 16 

2.2.3 Governing Equations 18 

2.2.4 Composite Laminate Element Stiffness Formulation ... 18 

2.2.5 Transient Analysis 21 

2.2.6 Explicit Integration Method 23 

2.3 Modeling of Damage 24 

2.3.1 Material Degradation 26 

iii 


PAGE BLANK NOT FILMED 



2.4 Analysis Procedure 27 

3 . Implementation 

3.1 OOP Concepts 28 

3.1.1 OOP Terminology 2 9 

3.1.2 Advantage of OOP 32 

3.2 EIFFE Class Libraries 34 

3.3 Class Descriptions 35 

3.3.1 CObject 35 

3.3.2 VN_ElementTemplate 36 

3.3.3 Communication between Element and Element 

Template Objects 37 

3.3.4 VN_Element 37 

3.3.5 VN_StructuralTemplate 3 8 

3.3.6 VN_3DTemplate 38 

3.3.7 VN_StructuralElement 39 

3.3.8 RKC_CompositeElement 39 

3.3.9 RKC_3DCompositeTemplate 40 

3.3.10 VN_FEModel 41 

3.3.11 VN_StructuralModel 41 

3.5.12 RKC_TransientModel 42 

4 . Case Studies 

4.1 Impact Response of a Steel Beam 44 

4.2 Transient Response of a Plate 48 

4.2.1 Isotropic Plate 48 

4.2.2 Composite Plate 51 

4.3 Inelastic Response of a Composite Plate 53 

4.4 Damage Prediction in a Composite Plate 56 


IV 



4.4.1 No Material Degradation 


56 


4.4.2 Material Degrdation 65 

4.5 Observations 71 

5. Conclusion and Recommendations 

5 . 1 Conclusions 72 

5.2 Recommendations 73 

References 74 

Appendix 79 


v 



LIST OF FIGURES 


Figure 

2.1 Impactor Striking Rectangular Plate at the Center 14 

2.2 Stacked Plies along the Thickness of the Element 19 

2 . 3 Location of Gauss Points for a 2 X 2 X 1 Integration 

Scheme in a Composite Ply 2 0 

3.1 The Top Level Class Hierarchy of EIFFE 3 5 

3.2 Class Hierarchy Element Template and Element Classes 36 

3.3 Node Numbering and Ply Numbering Scheme in a Composite 

Element 40 

3.4 The Class Hierarchy of Finite Element Model Classes 41 

3 . 5 Communication Process between Various Finite Element 

Objects 42 

4.1 Simply Supported Steel Beam 45 

4.2 Deflection at the Center of the Simply Supported Beam.... 46 

4.3 Finite Element Model of a Quarter Plate 48 

4.4 Transient Response at the Center of a Simply Supported 

Square Isotropic Plate subjected to Suddenly Applied 
Uniform Pulse Load 4 9 

4 . 5 Variation of Energies in a Square Isotropic Plate 

subjected to Suddenly Applied Uniform Pulse Load 50 

4.6 Deflection at the Center of a Square Composite Plate 

with a Lay-up Sequence [-45/45] subjected to Suddenly 
Applied Pressure Loading 51 

4 . 7 Deflection at Center of a Square Composite Plate with 

a Lay-up Sequence [30/45/90/0] subjecte to Suddenly 
Applied Uniform Pulse Load 52 

4.8 Deflection Response at the Center of a Clamped Composite 

Plate due to Inelastic Impact 54 

4 . 9 Variation of Energies in a Clamped Composite Plate due 

to Inelastic Impact 55 

4.10 Quarter Plate Finite Element Model of Clamped Composite 

Plate with lay-up Sequence [0 4 /-45 4 /45 4 /90 4 /45 4 /-45 4 /0 4 ] . . 58 


vi 



4.11 Deflection Response at the Center of the Composite 

Plate Model with No Material Degradation 59 

4.12 Resisting Force Profile at the Point of Impact in 

Composite Plate Model with No Material Degradation 60 

4.13 Variation of Energies in a Composite Plate Model 

with No Material Degradation 61 

4 . 14 Isometric View of the Predicted Damage with Model 

using No Material Degradation Method 62 

4 . 15 Top View of the Predicted Damage with Model using No 

Material Degradation Method 64 

4.16 Deflection Response at the Center of a Composite Plate 

Model with Material Degradation 65 

4.17 Resisting Force at the Point of Impact in Composite 

Plate Model with Material Degradation 66 

4.18 Variation of Energies in a Composite Plate Model 

with Material Degradation 67 

4.19 Isometric View of the Predicted Damage with Model using 

Material Degradation Method 68 

4.20 Top View of the Predicted Damage with Model using 

Material Degradation Mehtod 70 


vii 



TABLE 


TABLES 


4.1 Material Properties of Fiberite/T300/976 Graphite /Epoxy . . 57 


viii 



ABSTRACT 


Damage Tolerance 
in 

Filament- Wound Graphite/Epoxy Pressure Vessels 

by 

William E. Simon, Ph.D., P.E. Principal Investigator 
Vinh D. Nguyen, Ph.D., Research Associate 
Ravi K. Chenna, Research Assistant 

Graphite/epoxy composites are extensively used in the aerospace and sporting goods 
industries due to their superior engineering properties compared to those of metals. 

However, graphite/epoxy is extremely susceptible to impact damage which can cause 
considerable and sometimes undetected reduction in strength. An inelastic impact model was 
developed to predict damage due to low-velocity impact. A transient dynamic finite element 
formulation was used in conjunction with the 3D Tsai-Wu failure criterion to determine and 
incorporate failure in the material during impact. Material degradation can be adjusted from 
no degradation to partial degradation to full degradation. The developed software is based on 
an object-oriented implementation framework called Extensible Implementation Framework 


for Finite Elements (EIFFE). 
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CHAPTER 1 
Introduction 

1 . l Low-Velocity Impact Damage in Composites 

Graphite/epoxy composites are used extensively in the aerospace 
and sports equipment industries because of their high modulus and 
strength- to-weight ratios, extended fatigue life, and excellent 
corrosion resistance. These composites are replacing metals in 
aerospace structures such as aircraft wings and fuselages as well as in 
sports equipment such as bicycle frames, tennis rackets, and vaulting 
poles. Graphite/epoxy has been identified as a candidate material for 
pressure vessels in space applications (Lloyd and Knight, 1986) . 

However, graphite/epoxy is extremely susceptible to impact damage, 
especially at low impact velocities where internal damage may occur 
without manifestation on the surface. This particular kind of damage 
causes considerable reduction in strength and stiffness of the structure 
(Chang & Choi, 1991; Husman & Whitney, 1975; Jih & Sun, 1993; Poe & 
Garber, 1987) . For instance, experimental work done by Poe and Garber, 
(1987) on graphite -epoxy laminate cut from a filament -wound composite 
(FWC) pressure vessel showed 39% reduction in the laminate strength due 
to low-velocity impact damage. 

Several researchers have studied damage behavior of composite 
laminates undergoing low-velocity impact in the past two decades . These 
studies examined diverse topics including damage initiation, damage 
propagation , type of damage, influence of impactor mass and velocity, 
and laminate lay-up sequence (Kook et al . 1992; Preston & Cook, 1975; 
Chaturvedi & Sierakowski, 1985) . The studies on the effect of stacking 
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sequence on impact resistance indicated that laminates with more 
uniformly dispersed ply orientation have greater resistance (Choi et al . 
1992) . It was also found that impact damage is more sensitive to 
stacking sequence than to thickness. The damage zones were examined 
either visually, under an electron microscope, or acoustically using 
ultrasonic C-scanning and imaging (Chaturvedi & Seirakowski, 1985; Chang 
& Ketan, 1983) . The overall finding was that delamination accompanied 
by matrix cracking was the major damage mode associated with low- 
velocity impact. Delamination tends to occur at ply interfaces where 
the fiber orientation changes (Wu and Springer, 1988) . Choi and 
Chang (1991) , found that the delamination area is positively correlated 
to the impactor energy and recognized that out-of-plane tensile stress 
was the cause for delamination growth. Matrix cracking, on the other 
hand, usually occurs at the bottom-most layers. These matrix cracks 
were also found to initiate delamination at ply interfaces where the 
fiber orientation changed direction. 

1.2 . Experimental Studies 

Various experimental models including drop-weights (Madan, 1991; 
Oplinger & Slepetz, 1975), pendulums (Chou, Carleone & Mortimer, 1975; 
Cook & Preston, 1975), and air guns (Wu & Springer, 1986; Choi & 

Chang, 1991) . Impactors of different shapes and masses were designed to 
measure the contact force and impact velocities to determine the 
laminate impact response. Experimental results from both dynamic and 
static indentation tests suggested that under low-velocity impact 
conditions, the dynamic impact responses were similar to static 
indentation responses because the plots obtained for energy absorbed- to- 
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energy available ratio were similar to the force- indentation plots 
(Sjoblom et al . 1988) . Similar results were also obtained by Kwon and 
Bhavani (1993) when heavier impactor masses (1 to 10 Kg) with low- 
velocities (0-3 m/sec) were used. 

1.3 . Analytical Studies 

Simultaneously, analytical models were also being developed to 
simulate the impact behavior of the composite laminates. These models 
can be classified as either wave front models or contact models. 

1.3.1 Wave Front Models . Wave Front models simulate the impact 
load using approximation functions and Fourier series to describe the 
displacement field in the laminate. The equations of motion are derived 
in terms of these Fourier series and decoupled by transformation to 
modal coordinates (principal coordinate system) . The stress wave 
patterns, wave surfaces, and wave velocities in the material due to 
impact can then be obtained by solving the model equations (Kubo & 
Nelson, 1975; Moon, 1972, 1973; Slepetz & Oplinger, 1975). In these 
analyses, the loading behavior was simulated using certain approximation 
functions, whereas in reality, the impact loading behavior and the 
contact behavior is complicated and cannot be so easily represented. 

1.3.2 Contact Models . Contact models use Hertz ' s contact law as 
the basis for the computation of contact forces between the impactor and 
the impacted structure: 

F= na^/ 2 (1.1) 

where F is the contact force, 

a is distance between the centers of gravity of the two 


bodies, and 
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„ = 1 33*V^1Z^ + 1^£|, (1.2) 

^ E Ep ) 

Ep = Young's modulus of the projectile, 

E = Young 1 s modulus of the target , 

Vp = Poisson's ratio of the projectile, 
v = Poisson’s ratio of the target, and 
R = projectile radius. 

Hertz contact law predicted good results for isotropic materials, 
but failed for other types of materials. For transversely isotropic 
targets, Willis (1966) proposed using the transverse modulus of 
elasticity of the target in equation 1.2. Preston and Cook (1975), Sun 
(1977) used this approach to compute the impact response of a beam. 

These researchers solved non-linear integral equations to determine both 
the contact force and the impact response. While these studies 
demonstrated the possible application of Hertz's law to impact analysis 
of composite structures, the results obtained were not close to 
experimental results, probably due to the assumption that the impact was 
elastic . 

The modified Hertz law could not be successfully implemented for 
the following reasons: 

i) Most laminates could not be represented adequately by a 
half -space, 

ii) It did not account for anisotropic and nonhomogeneous 
laminate properties. 

Furthermore, while Hertz's law assumes perfect elasticity in both 
the impactor and the target, it was found that permanent indentations 
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may result at the point of impact. To account for permanent 
indentations, Crook (1952) , proposed the following equation for contact 
force : 


F = Fm 


a -a o 

OLm—CLo 


(1.3) 


where, F m = maximum contact force at unloading, 
am = indentation at unloading, 
ao = permanent indentation, 
q = unloading power. 

Yang and Sun (1982) found that indentation exist at very low- 
velocity impacts and that the loading and unloading behavior changes 
when maximum indentation exceeds a value called "critical indentation” . 
Contact parameters in the Crook's indentation law were derived by 
fitting experimental data was fitted for the loading and unloading 
processes. These indentation laws were successfully used by many 
researchers to determine the impact response and stress -strain histories 
(Kook et al. 1992; Chen & Sun, 1985; Sun & Tan, 1985). Tan and Sun 
(1985) successfully used them in a finite element model to study the 
dynamic response of a thin graphite /epoxy laminated plate with free -end 
boundary conditions and small deflections. The results obtained from 
this analysis agreed with experimental results. Sun & Chen (1985) 
investigated the influence of effects of inplane prestress on the plate, 
and the impactor's velocity, mass, and size. Previous studies were 
based on small -deflection theory which underestimated the magnitude of 
shear deformation. Later studies revealed that laminates undergo large 


deflections as well as transverse shear deformation during low-velocity 
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impact. Shivakumar et al . (1985) and Kook et al . (1992) included the 

effect of transverse shear deformation in their models. Shivakumar used 
thick plates (1/w < 12) in his study, while Kook used higher order shear 
deformation theory. Results obtained from these studies were close to 
experimental results. 

1 . 4 Residual Strength Prediction Models 

The most popular models available for predicting residual strength 
in composite structures are empirical. Such empirical models use 
fracture mechanics concepts to estimate the residual strength by 
establishing the relationship between residual strength, and impact 
parameters, such as the impactor velocity, or energy. The strength of 
any structure depends on its material integrity, i.e., if there is a 
flaw in the material, the strength of the material decreases. In such 
models, experimental data for the residual strength and the damage size 
is first compiled. An empirical relationships is then established 
between the governing parameters of damage, such as, the impactor 
velocity or impactor energy and flaw size. 

The residual strength of the structure is estimated in conjunction 
with experimental data and the established relationships. Caprino (1984) 
used a linear elastic fracture mechanics model to estimate the residual 
strength of composite laminates as a function of impactor energy. 
According to fracture mechanics principles, there exists a notch of 
characteristic length beyond which the structure experiences strength 
degradation according to the following formula: 



a 


(1.4 ) 
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where, a r = residual strength of structure for a notch of length C, 
a Q = residual strength of the structure for a 
characteristic notch length C 0 , 
n is the empirical constant. 

The damage caused is directly proportional to impact energy, i.e., 
the higher the impact energy, the larger the damage. This damage may be 
represented by a notch of an equivalent length determined as one created 
by the same amount of energy: 

C = kU™ (1.5) 

where, C = the equivalent notch length, 
k = Constant of proportionality, 

U = Energy required to create the notch, 
m = Empirical constant. 

Using the above relationship, the impactor energy for the 
characteristic notch can be determined. This equation, when substituted 
in equation 1.4, yields a relationship between residual strength and 
impactor energy: 



This model is simple and works well for materials which are rate- 
insensitive . 

Morton and Cantwell (1990), and Husman et al . (1975) represented 

damage as a sharp crack of equivalent transverse dimension. The 
residual strength was determined from experimentally established crack 
length- strength relationships. Poe and Garber (1987) determined the 
damage depth from analytical methods and represented this damage as a 
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notch, whose depth was determined empirically to determine the residual 
strength. The drawback to these methods was that damage was represented 
as a nonexistent hole or a crack and may not adequately describe the 
material degradation incurred. 

Tian and Swanson (1992) presented an analytical model consisting 
of damage analysis on a ply-by-ply basis using a 3D finite element 
model. This model accounted for fiber breakage and delamination as two 
separate failure modes. Line cracks were modeled to represent the fiber 
breakage and delaminations. A strain criterion was selected to predict 
the residual strength. The results obtained from this analysis closely 
approximated experimental data. 

1 . 5 Damage Characterization 

Damage characterization plays a central role in the prediction of 
residual strength as well as remaining life of composite structures. In 
the 1980s, the major thrust of research was in the area of damage 
quantification. Gu and Sun (1983) correlated specimen thickness change 
and reduction in the specimen's tensile modulus and strength. Ketan et 
al . (1983) concentrating on the point impact damage to sheet molding 
compound (SMC) panels, showed that the damage area is the primary 
parameter in characterizing the damage, rather than the indentation 
depth, or reduction of thickness. Gu and Sun (1987) proposed a 
semiempirical approach to predict the damage zone in SMC panels using 
the strain energy density function. The method is based on the 
assumption that composite failure is caused by excessive strain energy. 

A weighted strain energy function was proposed: 
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U = U v + pU s (1.7a) 

where U v = dilatational energy, 

U s = distortional energy, 

P = is a weighting coefficient. 




1 - 2v 


( CT XX xx® yy + ° yy) 


2(1 + v) 7 j 7 ~ 

Us= “ (g xx ~ & xx® + ® „vv) + 2(1 + v Xcr xy + cr yz + G xz) 


(1.7b) 

(1.7c) 


The strain energy was expressed in terms of equivalent stress a 
as described by Equation 1.8 during the analysis. 

a 2 = 2EU (1.8) 

where, E = modulus of elasticity, and 
U = strain energy. 

Failure was assumed to occur at a critical value of U, U cr or at 
critical value of a, a cr . A transient finite element analysis was 
conducted with Mindlin's plate elements. The effective stresses for the 
entire analysis was stored. The maximum value of a at the Gauss points 
yielded a distribution of effective stresses. To match the analytical 
solution with experimental data, two variables: the critical stress 

a cr and the weighting coefficient p, were chosen to fit the 
experimental data. The same values were used to predict the damaged 
areas for different cases of impact velocities. It was found that the 
value of P was usually less than 1, implying that dilatational energy 
(energy associated with dilatation strain) plays a dominant role in 
failure mechanism. This failure criterion proved to be adequate for a 


wide range of impact velocities. 
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It was shown in many other studies that matrix cracking and 
delamination were the major modes of damage in low- velocity impact 
(Chaturvedi & Sierakowski, 1985; Choi & Chang, 1991; Jih & Sun, 1993; 
Madan, 1991) . Wu and Springer (1986) revealed that out-of -plane stress 
is the major cause of delamination. These authors developed a three- 
dimensional transient finite element model of a simply supported plate 
subjected to low-velocity impact. The Tsai-Wu failure criterion was 
used to predict delamination damage which was shown to agree with 
experimental data. No material degradation was included in this 
analysis. Later Wu and Springer (1988) developed a model based on the 
dimensionless parameter theory relating the delamination dimensions with 
parameters influencing the size of damage or delamination. 

1 D = f (CT,R,tf ,AQ ( D t ( D b , 1 0 ,K c ) (1.9) 

where, a = is the stress at the location of the damage, 

R = is the rate at which the stresses change, 
t f = is the duration of the stresses, 

A Q = is the difference in reduced stiffness of the two 
plies adjoining the delamination, 

D T 'D b = are the flexural rigidities of the layers above 
and below the interface where delamination occurs, 

1 0 _ initial size of the flaw and 

K c = is the resistance of the material to separation. 

The predicted delamination dimensions were within 20% of the 
experimental results. Again, material degradation during the impact was 
not considered. This may be the reason for the difference between the 
analytical results and experimental data. 
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Choi et al. (1991) conducted experimental studies to show that 
matrix cracking and delamination play a major role in impact damage. 
They used a line-nosed impactor to produce a uniformly distributed 
transient dynamic load across the specimen* s width. The following 
conclusions are drawn: 

i) matrix damage was the initial damage due to impact; 

ii) delamination follows the "critical" matrix crack; 

iii) impact energy threshold is governed by the energy required 
to initiate the first critical matrix crack; 

iv) stacking sequence affects impact resistance of composites. 

The experimental studies by Choi et al . (1991), showed that 

interlaminar shear stresses and in-plane tensile stresses were dominant 
factors causing matrix cracking. These matrix cracks produce micro- 
cracks which in turn caused delamination failure. The out-of -plane, 
stresses on the other hand cause delamination growth. Choi et al . 

(1991) developed a two-dimensional transient finite element model to 
verify that matrix damage indeed initiates delamination. A modified 
Hertzian contact for line- loading was used in conjunction with a matrix 
failure criterion. The material stiffness of the elements within the 
damaged area was degraded and the stresses redistributed. Choi and Chang 
(1992) also developed a three-dimensional transient analysis to verify 
that out -of -plane stresses cause delamination growth. In this model, in 
addition to the matrix failure criterion, a delamination failure 
criterion was also included. To predict damage, the material was first 
tested for matrix failure and then for delamination, at the locations 
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where the matrix failed. This model was capable of predicting both 
matrix cracking and delamination due to the impact. 

1.6. Outline pf present work 

In the present work, another model is proposed to predict damage 
due to low-velocity impact. The proposed approach is similar to that 
used by Wu and Springer (1986) , except that degradation is included 
during the impact. Also, instead of the Hertz contact model, an 
inelastic impact was used which is simpler and effective. To predict 
the extent of damage, the model utilizes a generally accepted 3D Tsai-Wu 
failure criterion. While the previous studies did not include material 
degradation during impact the present model allows for material 
degradation at regular intervals . The model was implemented using 
object oriented programming (OOP) concepts to support future extensions. 
Consequently, the software can be easily modified and specialized for 
different conditions and applications (e.g., incoporating Hertz's 
contatct law) without much difficulty. 

1.6.1 Or ganization . Chapter 2 discusses the mathematical 
formulations for modeling impact as well as the composite element 
formulation. Chapter 3 discusses the structure of the software package 
developed to implement the model . The results and discussions are 
presented in chapter 4 . Conclusions and recommendations are presented 


in chapter 5 . 
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CHAPTER 2 
Formulations 

This chapter discusses the mathematical formulation and concepts 
used to develop the finite element model of impact damage in composite 
laminates. The approach adopted is similar to the one used by Wu and 
Springer (1986) , with one important exception. The Wu- Springer model 
did not include material degradation during damage, whereas the present 
work includes cumulative degradation which can be applied at user- 
specified intervals, or not applied at all. Section 2.1 provides a 
brief problem statement. Section 2.2 discusses the various formulations 
used in the model. Section 2.3 lists the step-by-step analysis 
procedure . 

2 . l Statement p£ the problem 

The objective of the present work is to develop a finite element 
model to predict damage in composite structures undergoing low velocity 
impact and accounting for cumulative material degradation during impact. 
The model is to be validated with a case study involving a rectangular 
plate whose dimensions are L (length) , W (width) , and h (thickness) , 
with continuous fibers. The plate consists of n plies whose 
orientations are arbitrary and need not be symmetric to the midsurface 
of the plate. An impactor of mass Mp strikes the stationary plate with 
a velocity V p as shown in Figure 2.1. While developing the finite 
element model following conditions were assumed: 

i) perfect bonding between plies; 

ii) individual layer are homogeneous and orthotropic, and 

iii) the impactor adheres to the target after impact. 
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Figure 2 . 1 

Impactor Striking the Rectangular Plate at the Center 
(Mp: Mass of the Impactor, Mn: Mass of Impacted Node) 


2.2 Analysis 

The equation of motion in variational form can be expressed, 
without damping, as: 

l n w, p u ltll dv + l n ejj Ejjue k/dv - l r WjOjjnjdA = 0 (2.1) 

where, = the stress tensor on the boundary of the domain, 

e kl = the strain tensor, 
p = the density of the plate material, 
u iitt = acceleration vector, 
w i = arbitrary variations, 

Q = the volume of the domain, 

T = the domain boundary. 


nj = the outward normal vector on the domain boundary, 
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E ijki = fourth order tensor represents the material 
stiffness . 

The three terms in Equation 2.1 represent the kinetic energy, the 
strain energy, and the work done by the boundary forces, in this case 
the contact forces . 

2.2.1 Modeling Contact Behavior. Hertz contact model is a 
frequently used to predict contact forces. As suggested by Lin and Lee 
(1989) , impact may be modeled as inelastic when the impactor mass is 
larger than the mass of the impacted node. The authors found that the 
results obtained from such a model were in good agreement with 
experimental data. As a result, impact may be modeled as an initial 
velocity at the impacted node with conservation of momentum taken into 
account . 


A similar inelastic contact model is used in this present work. 
The conservation of momentum principle dictates that the equivalent 
initial velocity of the impacted node is given by: 


Vo = 


MpVp 


where, Mp = the impactor mass, 


( 2 . 2 ) 


Vp = the impactor velocity, 

Mn = the mass of the impacted node on the plate, 


Vo = the equivalent initial velocity of impacted node. 
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2.2.2 Finite Element Formulation. An 8-noded brick element was 
used to model the composite plate in a transient finite element 
formulation. The displacement at any point in the laminate can be 
expressed as: (Wu & Chang, 1987) 

8 

Uq ~ ^ N r Uq r q = 1 , 2 , 3 # 2.3 

r-\ 

where, u qr are the nodal displacements, and N r (r = 1 ~ 8) is the 
shape function vector for a 8-noded brick element, N r can be written as 
follows : 

N, » (1-0 (1-T1) (1-0 / 8 , 

N 2 = (1+0 ( 1 -t| ) (1-0 / 8 , 

N 3 = (1-0 (1+Tl) (1-0 / 8 , 

N 4 = (1+0 (1+11) (1-0 / 8 / 

N 5 = (1-0 ( 1 -T 1 ) (1+0 / 8 , 
n 6 = (i+o (i--n) u+o / 8 , 

N 7 = (1-0 (1+Tl) (1+0 /8, 

N 0 = (1+0 (1+Tl) (1+0/8, (2.4) 

where 0 r| and are the natural coordinates of the element. 

An isoparametric formulation is used such that the coordinates 
(q = 1-3) of any point inside the element can be expressed in terms of 
the shape functions : 

8 

Xq = N r Xqr q=l,2,3 (2.5) 

r-\ 

where x qr are the coordinates of the node r. 

The strains at any point in the laminate can be expressed as 
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[e] T = Z[fi r ]{^| r U 2 r U ir f 
r= 1 


where 


( 2 . 6 ) 



N r _] 

0 

0 

0 

N r . 3 

N r . 2 

[Br] = 

0 

N r . 2 

0 

N r . 3 

0 

N r . 1 


. 0 

0 

Nr. 2 

N r . 2 

Nr.] 

0 


8 , 


The stresses at any point are given by the following relation 
CT ij = Qijki E ia i,j,k,l = 1-3 (2.7a) 

where Q ijkl is material stiffness tensor, which varies with fiber 
orientation. The above equation can be reduced to Equation 2.7b by 
applying material symmetry. Detailed explanation can be obtained from 
any text book on continuum mechanics or any text book dealing with 
mechanics of composite materials. 

tfi = Qij Ej i, j = 1-6 (2 . 7b) 

Qij for a lamina is computed along its principal material 
direction i.e., in the direction of fiber, and then Q # ^ is obtained by 
rotating the Q by the fiber orientation. The following Equation 2.7c 
shows the transformations required to rotate thru angle 0. 


[Q'] = [T] ’ 1 [Q] [T] 

where, [Q' ] = material stiffness matrix in the global coordinate 
system, and 


[T] = transformation matrix. 
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2.2.3 Governing Equations . The finite element equation for 
transient analysis is given by Bathe & Wilson, 1976. 

[M]{1}+[C] w+[/:] {*} = {*} ( 2 . 8 ) 

where, [M] = mass matrix, 

[K] = stiffness matrix, 

[C] = viscous damping matrix (a [M] + P [K] ) , 
a and p are experimentally determined constants, 

{ X } = nodal displacement vector, 

{R} = external load vector. 

The basic components involved in the above equation are the mass 
and stiffness matrices. For computational convenience, a diagonal 
lumped mass matrix was used. In this method the mass of the element is 
equally distributed at its nodes. The expression for mass matrix is: 

[M e ]=—[I] (2.9) 

8 

where, [I] is 24 x 24 identity matrix. 

2.2.4 Composite Laminate Element Stiffness Formulation. The 
formula for stiffness matrix of a finite element is given as 
(Zienkiewicz , 1977) : 

[K e ]=\ ne [B] T \Q][B)dV ( 2 . 10 ) 

where, [Q] is the material stiffness matrix. 

In composite laminates, [Q] varies from layer to layer due to the 
change in the fiber orientation. There are two methods available to 
model the composite element. In the first method, each ply is 



Chenna 19 


represented by an element along the thickness of the laminate. This 
approach results in a large model, making the analysis computationally 
expensive. In the second approach, each composite element is assumed to 
contain stacked plies oriented at different angles along the thickness 
of the element, as shown in Figure 2.2. The second approach allows the 
material properties to vary in discrete steps through the thickness of 
the element . 



Stacked Plies along the Thickness of the Element 

The stiffness of a composite element is computed by summing the 
stiffnesses offered by each ply along the thickness direction of the 
element : 



( 2 . 11 ) 
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where, N is the number of plies along the thickness of the 
element, Z k+1 and Z k are the top and bottom coordinates of the kth ply in 
the element. [Q] k is the material stiffness matrix of the kth ply. The 
global coordinates x, y, z are mapped into a natural coordinate system 
having coordinates ^,r| and C, {Zeinkewicz, 1977). 

The stiffness offered by each ply in an element is computed by 
the summation of the stiffnesses at the Gauss points in the ply. The 
Gauss points in a ply with 2X2X1 integration scheme is shown in 
Figure 2.3. 



Figure 2 . 3 


Location of Gauss Points for a 2 X 2 X 1 Integration 
Scheme in a Composite Ply. 
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2.2.5 Transient Analysis 

At any time t, the load vector { R } in Equation 2.8 can be written 
as a sum of the following components; 

Fi( 0+ Fd(0 + FeU) = F(t) (2.12) 

where, F j{t) = [M]{X} t are the inertial forces, 

Fo( t )-[C]{X} t are the damping forces, 

Fe(0 = [K]{X} t are the elastic forces. 

All of the above variables are t ime- dependent . Mathematically, 
Equation 2.8 is a second order differential equation which is solved 
using direct integration. In the direct integration method, the 
derivatives are approximated by finite differences and integrations may 
be carried out implicitly or explicitly. 

In the implicit method, the equation of motion is satisfied at the 
next time step, i.e., at t+At. The following equation illustrates the 
Newmark method: 

[mX} t+ ^[C]{X} t+ ^[K]{X) l+ ^ = {*}, +A/ (2.13) 

W, +A/ = W,+(1-y)^W,+yA/{1} /+A/ (2.14) 

W, +A/ =w, + a/w, + (03-P)a/ 2 w,+ p a/ 2 w /+ a , (2.i5) 

where, (3 and y are the experimentally determined parameters. 

Using the above Newmark finite difference Equations 2.14, and 2.15 

for displacement and velocity vectors at t+At, the terms {X} (+ ^ t and 

{X} /+ a / can be derived in terms of {X} t+diI . These expressions, when 
substituted in the Equation 2.13 above, result in the following equation 
expressed in terms of displacement vector {X} t+Al at time t+At : 
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[^]W,+Af = W, +A , (2.16a) 

where, [AT] is the effective stiffness matrix defined as: 

K = K + a 0 M + a\ C (2.16b) 

1 8 

a 0 = ; a\ = 5 2: 0.5; a ^ 0.25(0.5+8) , 

a At 2 a A t 

and 



5 AffS 

#4 = — - 1 ; a 5 = — — -2 

a 2 \a 

The solution of Equation 2.16a is substituted back in the finite 
difference equations to compute the velocity and acceleration vectors at 
time t+At. There are many other finite difference methods available, 
namely the Houbolt method, the Wilson 0 method, and Central difference 
method. The procedure involved in all these methods is similar (Bathe & 
Wilson, 1976) . The advantage of implicit method is that large time 
steps can be used. The disadvantages of implicit methods are 
(Belytschko, 1983) : 

i) The methods are computationally expensive since they involve 

the assemblage and inversion of [A*] . Also they are not 
suitable for structural problems involving progressive 
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failure because of the numerical difficulties that often 
result from material degradation (Wolcott & Yener, 1989) ; 

ii) They demand greater storage space because of the complexity 
and size of software; 

iii) They cannot be applied easily to nonlinear problems as the 
changes that must be made to the software make the methods 
even more computationally expensive. 

2.2.6 Explicit Integration Method. In this method the equation 
of motion is satisfied at the current time step t, i.e., 

[A'lWf + ICKAVMW, -<*}/ <2 . 17a) 

{X}, = [ M]~ l {{*}, -[C]{X) t ~[K]{X} t (2 . 17b) 

As inelastic impact and no damping were assumed the equation 2.17b 
may be expressed as: 

{X}, =[A/r'{-[K]W, (2 . 17c) 

Using the present acceleration vector, the velocity and 
displacement vectors can be derived for the t+At time step using simple 
finite difference equations: 


W, +Af -{*},+{*},* A/ 

(2.18) 

{X} t+ ^{X} t + {X)*M 

(2.19) 


The displacement and velocity vectors computed at time step t+At 
are used to update internal resisting forces. This procedure is carried 
out for the entire impact time duration to obtain the complete response 


of the system. 
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Following are the advantages and disadvantages of the explicit 
integration methods (Belytschko 1983) : 

i) Fewer computations are required in each time step; 

ii) The algorithm is simple in logic and structure; 

iii) They requires little core space compared to implicit 
methods ; 

iv) They are ideal for testing new ideas, because it requires 
less coding; 

v) They are conditionally stable, so that very small time steps 
may be required. 

The small time steps are compensated for by simple computation 
required in each time step. To assure stability, the time step must 
always be less than a critical time step computed from the largest 
eigen- value of the system. Critical time is computed using (Belytschko, 
1983) : 



where At cr is critical time step and 'k TMX is the maximum eigen 
value computed on an element -by- element basis. The explicit method was 
adopted for present work. 
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2.3. Mpdeling._Qf Damage. 

Failure in composites is complicated by a multitude of interacting 
mechanisms including fiber breakage, micro-buckling, delamination, fiber' 
pullouts, and matrix cracking. The Tsai-Wu failure criterion (Tsai, 

1971) is used to determine overall damage without distinguishing 
different types of failure. The failure criterion for orthotropic 
material is expressed by the following inequality: 

F a a x + F 2 a y + F 3 a z + F 4 a y2 + F 5 a 2X + F^ 

+ Fn<*x l 2 + 2F 12 CT x a y + 2F 13 o x a z + F 22 a y 2 + 2F 23 a y a z + F 33 a z 2 
+ F 44 a yz 2 + F 55 a zx 2 + F^a^ 2 ^ 1 (2.21) 

where, ct x , a y , a 2 are stresses along X, Y, Z axes, respectively and 


°yz' CT zx/ a xyr are shear stress in Y-Z, Z-X and X-Y planes, 
respectively. 

The coefficients F 1# F tj are given as: 

1 1 1 

F\ = — - — , Fn = , 

X X * XX’ 


1 1 1 

F 2 = “"“/ F 22 = 

Y Y YY 


1 1 


2 Z 

1 1 

F 4 = — — 

Q Q 


l 


F 3 = ~ / ^33 = ~f 


F 44 = ' 


zz 

\ 

QQ 


l l 

f 5 = — — 

R R’ 


^55 = —' 
RR 


1 1 1 

F 6 = — / F$6= / 

" ~ SS' 


F \2 = ~ 
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where, X and X' are tensile and compressive strengths along the 
fiber or X, direction, 

Y and Y' are tensile and compressive strengths of the 
material along transverse, or Y, direction, 

Z and Z ' are tensile and compressive strengths of the 
material along normal, or Z, direction, 

Q and Q' are positive and negative pure shear strengths on 
the Y-Z plane; R and R' , on the Z-X plane; and S and S' 
along X-Y plane, and 

P is experimentally determined by application of biaxial 
tension. 

2.3.1 Material Degradation. Material degradation can be either 
be complete or partial. In the complete degradation method, the 
stiffness at failed Gauss point is completely removed from the overall 
element stiffness: 

[K'] = [K] prev - [K] failcd point8 (2.22) 

where [K' ] is damaged stiffness matrix, 

[K] prev is the previous stiffness matrix, 

[K] failed points the stiffness at the failed points. 

In the partial degradation method, only a fraction of the 
stiffness at the failed Gauss points is removed: 
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IK'] = [Klprev - <*>* [K] failed points (2.23) 

where O is the degradation factor (o £ £ 1) . 

2 . 4 Analysis Procedure 

The following steps summarize the analysis procedure: 

i) A 3-D finite element model is developed. 

ii) Transient analysis of the finite element model is carried 
out to determine the impact response. 

iii) During the analysis, nodal displacements and internal 
resisting forces are computed for every time step At. 

These displacements are used to compute stresses and strains 
at the Gauss points. Also, strain energy, kinetic energy, 
and total energy are computed, and recorded for the whole 
model . 

iv) The computed stresses at the Gauss points are checked with 
the Tsai-Wu criterion for material integrity. 

v) The Gauss points at which the material fails are recorded, 
i.e., the Gauss point coordinates. Depending on the 
material degradation method, the element stiffness may be 
either completely or partially degraded. 

vi) After material degradation, the acceleration vector at 
current time step is computed. Using the accelerations, the 
displacements, velocities and the internal resisting forces 
at step t+At are determined. 

vii) Steps iii thru vi are repeated for the duration of interest, 

viii) The failed Gauss point locations recorded during the 


analysis are used to plot the damage regions. 
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CHAPTER 3 
Implementation 

The software package developed is based on an object-oriented- 
implementation framework called Extensible Implementation Framework for 
Finite Elements (EIFFE) . Section 3.1 discusses the fundamentals of 
Object Oriented Programming (OOP) concepts. Sections 3.2 and 3.3 
describe the various EIFFE classes as well as classes developed to 
implement the model described in chapter 2 . 

3.1 qqp Con c epts 

Early engineering software was written in procedural languages 
that modeled logical units as black boxes. Every unit of code is boxed 
off so that it remains concealed. Such boxes are called functions in C 
and procedures in Pascal (Eckel, 1993) . The major drawback of 
procedural programming is that it lacks control over data. The 
designers of procedural programming languages designed the languages 
based on the assumption that the code required no maintenance or 
extension. Such assumptions worked when the project size was small, but 
failed miserably when the complexity of the problem increased. 

To accommodate the complexity of large programs, the concept of 
structured programming was introduced. In structured programming, large 
programs are divided into modules. Each module is further divided into 
sets of related procedures that manipulate data. Although structured 
programming improved clarity, reliability, and ease of maintenance of 
software, large-scale programs continues to pose challenges, due to the 


following reasons: 
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i) Structured programming require a great deal of planning to 
generate software that is extensible, maintainable and bug- 
free . 

ii) The software developed is usually rigid and intractable; 
modification is difficult due to the interaction of code 
from different modules. 

In, the 1980s, a new programming paradigm called OOP evolved, 
which alleviated some problems. While procedural programming hides the 
complexity of operations performed on data, OOP hides both the data and 
the operations {Eckel, 1993). An object-oriented language emphasizes 
data types and the intrinsic operations that may be performed on those 
data types . Data do not flow openly around a system as in procedural 
programs, but are protected from accidental modification. In OOP, 
function calls are replaced by message passing. Messages cause objects 
to manipulate data. 

3.1.1 OOP Terminology. The following are some terms often used 
in object oriented languages. 

Object : An Object is an instance of a class, is defined by a set 

of attributes, (e.g., a geometrical object would have points, color, and 
size as the variables) and a set of procedures (rotate, move, expand, 
contract) that operate on these attributes. 

Class : A class is a template for creating objects which share 

common attributes, and methods. A class not only defines the object's 
attributes but also provides methods for manipulating these attributes. 
The following example defines a class called CPerson: 
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class CPerson 

{ 

private : 

CDate m_Da t eof Birth ; 
float m_Age,- 

char* m_Name / m_SSNumber, m_Address, m_Telephoneno; 
public : 

void Person (char* Name, char* SSNumber, CDate DateofBirth) ; 

void EnterAddress (char* Address); 

void EnterTelephoneNumber (char* TelephoneNo) ; 

float GetAge ( ) ; 

char * GetAddress ( ) ; 

char * GetTelephoneNo ( ) ; 

virtual void DisplayData ( ) ; 

virtual void PrintDataO; 

}; 

class CDate 

{ 

private : 

int month, day, year; 
public : 

SetDate(int, int, int); 


GetDate(int *, int *, int *) ; 
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Data Encapsulation: The isolation of data from the external 

environment is called data encapsulation. This feature is also known as 
data hiding. Since the data stored in the object is inaccessible, their 
modification can be performed only through a controlled interface. Data 
Encapsulation promotes modularity. In the above example, data such as 
name, social security number, address in class CPerson cannot be 
accessed directly. However, they can be manipulated using CPerson' s 
interface, i.e., by calling its methods. 

lube r i t..anc.s i Inheritance defines a relationship among classes, 
wherein a class inherits the attributes and behavior of one or more 
other classes. This feature allows the software to be reusable. The 
following class CStudent example illustrates the principle of 
inheritance and the software reusability feature of OOP. 
class CStudent : public CPerson 
{ 

private : 

char* m_University , *m_ClassNo; 
float m_GPA; 
public : 

CStudent (char* Name, char* SSNumber , char* Dtof Birth, char *Univ, 

char* ClassNo) : CPerson (Name , SSNumber, DtofBirth) { } ; 

char * GetUniversity { ) ; 

char * GetClassNo ( ) ; 

float GetGPAO; 

void DisplayData (COutput OutputObject) ; 

}; 
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Attributes such as name, social security number, and date of birth 
are inherited from class CPerson. The relationship can be clearly 
understood if student is seen as a kind of person with some additional 
attributes. Other attributes, such as the name of university, and 
classes the student is attending, are added to the CStudent class. As 
most of the attributes and methods are inherited, little work is needed 
to create the class CStudent. 

Polymorphism; It is a concept wherein different objects may 
respond differently to the same message depending on the classes the 
objects belong to. A message "DisplayData" in the above example, sent 
to a computer console object results in data printed on the computer 
monitor, whereas the same message sent to printer object results in data 
being printed on paper. The message is the same but the response is 
different . 

3.1.2 Advantages of OOP. OOP provides better security and 
reusability than procedural programming in many ways: 

i) Objects are self-contained entities that can be introduced 
into a system without much difficulty, since any bug 
associated with the new code is localized to the object 
itself . 

ii) New object types may be derived from previously defined 
ones. This saves time and supports quick explorations. New 
objects do not modify the behavior of parent objects and 
keep new bugs localized in the new objects. 

iii) Well designed classes support code reusability and 


extensibility and lead to greater productivity. 
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iv) Partitioning of work comes naturally, thus making delegation 
of work in large projects easier. 

v) Data encapsulation promotes secure systems. 

vii) Software maintenance and management becomes much easier due 
to encapsulation and uniform object interfaces. 

viii) Data Encapsulation increases readability and reduces the 


need for documentation. 
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3 .2 EIFFE Class Libraries 

EIFFE consists of five class libraries: 

i) Matrix classes, 

ii) Function classes, 

iii) Finite element classes, 

iv) Material classes, 

v) Finite element collection classes. 

The first level organization is shown in Figure 3.1. Various types of 
matrix objects such as symmetric, banded symmetric, rectangular 
matrices, and vectors may be created using the matrix class library. 
Basic matrix operations such as matrix addition, subtraction, 
multiplication, and several other specialized techniques for solving 
matrix equations are also defined. Different types of functions such as 
monomial, bivariate, trivariate, and polynomial function can be created 
and evaluated at any point with the function objects created from the 
function class library. The finite element class library defines 
various finite element objects such as nodes, boundary conditions, and 
elements that are required to build a finite element model. The type of 
material the model uses can be defined as anisotropic, orthotropic, 
transversely isotropic, or isotropic material objects defined in 
Material class library of EIFFE. The methods for computing different 
kinds of stiffness matrices such as plane-stress, and plane-strain 
stiffness matrices are also defined. Material failure computations using 
the Tsai-Hill criterion and Tsai-Wu criteria are included. Finite 
element objects can be stored in collection objects, such as element 


list, node list, and function list. The collection classes are included 
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in the finite element collection class library. A brief description of 
finite element classes used to develop the present software is 
described . 




Figure 3 . 1 


The Top Level Class Hierarchy of EIFFE 
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3.3.2. VN_ElementTemplate The VN_Element Template class has 
methods defined to compute shape functions, shape function gradients, 
and the Jacobian matrix. These methods are used in the computation of 
strain-displacement matrices, stiffness matrices, internal or initial 
force vectors, and element stress and strain vectors. Since different 
element types may perform some of these computations differently, these 
methods may need to be implemented or overridden in the derived classes 
that represent particular finite element formulations. The template 
object is invoked by an element object to perform certain computational 
tasks (computation of element stiffness matrix, initial force vector, 
etc) . Before other types of element template classes are explained, it 
is important to understand the relationship and communication process 
between an element object and element template object. 



Figure 3 . 2 


Class Hierarchy of Element Template and Element classes 






Chenna 3 7 


3.3.3 Communication between Element and Element-Template objects 

When an element object is asked to compute its matrices, i.e., 
stiffness matrix, initial force vector or mass vector, it forwards the 
message to the attached element template object. The template object, 
before processing any request, checks whether the requesting element is 
of the right kind. The computational methods vary with the type of 
element (2D or 3D) , and the element order, that is, the order of shape 
functions (linear, quadratic, etc) . The template object also checks 
whether the requesting element has the right number of nodes to ensure 
compatibility between shape functions specified for the template and the 
number of nodes in the element, e.g., Four nodes are required for 2D 
element with linear shape functions, nine nodes are required if shape 
functions are quadratic. The template object uses the nodal information 
provided by the element object to compute the Jacobian matrix. It then 
requests the material object to compute the material stiffness matrix 
which is required in computation of the element stiffness matrix and the 
initial force vector. Once the element template object has all the 
required components, it creates the appropriate function integrands and 
requests the suitable integrators to carry out integrations . The 
computed matrix or vector is finally returned to the requesting element. 
Figure 3.2 shows the communication process between various element 
objects and element template objects. 

3.3.4. VN_Element ; The VN_Element class defines element objects. 

It includes the basic methods for computing stiffness matrices and 
stress and strain vectors. These methods simply call the corresponding 
methods in the template classes. Each element object also encapsulates 
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the connectivity data needed to evaluate the Jacobian matrix. The 
element object assembles its matrices, and vectors, into the global 
matrices and vectors . Similarly, it retrieves element matrices and 
vectors from global matrices and vectors. The hierarchies of element 
classes and element template classes are shown in Figure 3.2. 

3.3.5. VN_St rue tural Template : VN_Structur a 1 Template is a 

subclass of VN_ElementTemplate that deals with structural analysis. The 
methods defined deal with computation of matrix and vector objects, such 
as, stiffness matrix, body force vector, initial force vector, mass 
matrix, stress and strain vectors, etc. The structural element template 
object uses suitable integrators to integrate different types of 
function integrand objects provided to it. Using different element 
formulations, function integrand objects are created by VN_3DTemplate , 
VN_PlaneStrainTemplate, and VN_PlaneStressTemplate, which are subclasses 
of VN_StructuralTemplate class . 

3.3.6. VN_3DTemplate : This class is a subclass of 

VN_Structural Template . All the computational methods required by a 
3-dimensional element are defined in this class. This class is capable 
to compute the strain-displacement matrix [B] , for a 3 -dimensional 
elements. It uses the [B] matrix and material stiffness matrix, [Q] , to 
create the suitable 3D integrand functions such as stiffness integrand 
function, initial force integrand function, etc. The order of the 
functions generated depends upon the shape function generator used when 


creating the template object. 
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3.3.7. VN_StructuralElement : The VN_StructuralElement class 

inherits from the VN_Element class. VN_StructuralElement objects are 
used in finite element models dealing with structural problems. 
Structural problems can be solved using any one of the finite element 
formulations such as plane strain element, plane stress element, or 3D 
brick element. A structural element object can thus use any of the 
template objects created from VN_PlaneStressTemplate, 
VN_PlaneStrainTemplate or VN_3DTemplate . Structural element objects 
delegate the task of computing their matrices and vectors to the 
appropriate structural template object. The template object performs 
the computations and returns the results to the structural element 
object . 

3.3.8. RKC_CoTtpositeElement : The RKC_CompositeElement inherits 

from VN_StructuralElement . A composite element object is assumed to 
consist of multiple plies stacked along the thickness direction. The 
composite element object keeps track of the thickness and the lay-up 
sequences of the plies. This information is passed to the composite 
template object upon request. The node and ply numbering scheme used, 
are shown in Figure 3.3. The composite element also tracks the damage 
status at the Gauss points. When the composite element is asked to 
update its stiffness, it in turn requests the composite template object 
to compute the stiffness degradation at the damaged points. The element 
object either completely or partially removes the stiffness at the 
damaged points depending on the degradation method chosen. 
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c 


Figure 3 . 3 



Node Numbering and Ply Numbering Scheme in a Composite Element 


3.3.9. RKC_3DCompositeTemplate : This class is similar to the 

VN_3DTemplate class, with difference being that this class uses an 
orthotropic or transversely isotropic material object to model the 
composite material properties. The composite template object requests 
the material object to computes its material stiffness in the fiber 
coordinate system and then transforms it to the global coordinate 
system. It also has methods to compute stresses and strains in various 
plies . Methods are defined in this class to check for failure in the 
material using the Tsai-Wu criterion. After checking for failure it 
updates the material status in the element object by keep track of 


failed points . 
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Q VN_FEModel ) 


^VN 

_StructuralModel ^ 



^ RKC_TransientModel ) 


Figure 3 . 4 

The Class Hierarchy of Finite Element Model Classes 

3.3.10. VN_FEModel : VN_FEModel is the manager of the whole finite 
element model. VN_FEModel stores and keeps track of all the elements, 
nodes, and boundary condition objects belonging to the finite element 
model. Finite Element related objects can be added or removed. 
VN_FEModel has a GO method that instructs the model to proceed with the 
analysis. The hierarchy of finite element model classes is shown in 
Figure 3.4. 

3.3.11. YN__£ true t uralModel ; This class is a subclass of 
VN_FEModel that handles structural analysis models. The VN_Structural 
Model object processes user request by sending its own message to the 
objects it contains in the correct order. For example, after the 
creation of a finite element mesh, each element object in the model is 
requested by the structural model object to its own element stiffness 
matrix. When the structural model is requested to solve the model, it 
sends messages to element objects to request the assembly of their 
element matrices and vectors. After assembly, the model object issues a 
message to the appropriate solver object to obtain the analysis results. 
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3.3.12. RKC_Transien tModel : This class inherits from 

VN_StructuralModel . The transient model object is used for solving 
transient problems wherein the solution is sought at different time 
marks. The methods involved could be classified into four stages, the 
different stages are shown in Figure 3.5. 

Stage 1 : The transient model object calls all the elements to 

compute their degrees of freedom, mass vectors, and stiffness matrices. 
The element objects in turn requests its element template object to 
compute the element vectors and matrices. The template object send 
message to material object to compute the material stiffness matrix 
which is used in the computation of stiffness matrix. The computed 
element matrices and vectors are stored in element object. 

Stage 2 ; The transient model object creates empty vectors for 
nodal displacements, velocities, accelerations and internal resisting 
forces. It also sends message to each element object to assemble its 
element vectors into global vectors. Both the initial conditions and 
the essential boundary conditions are also applied. 

Stage 3 ; The equations of motion are solved for nodal 
accelerations . 

Stage 4 : The velocity and displacement vectors are updated for 

the next time step. The elements are requested to update their state, 
i.e., check and record any material failure as well as updating the 
stiffness matrix. Internal resisting forces are recomputed using 
updated stiffness matrix and displacement vector. Stages 3 and 4 are 


repeated for the entire duration. 
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CHAPTER 4 
Case Studies 

Four case studies were used to validate the model and to 
illustrate the analytical capabilities built into the model. In 
section 4.1, the impact response of a simply supported steel beam was 
studied. The composite element formulation is verified in section 4.2 
by solving the transient response of a steel plate and a composite 
plate, both subjected to a suddenly applied uniform pulse load. The 
i n ®isstic impact response of a composite plate is verified in 
section 4 . 3 by solving a numerical example taken from Lin and 
Lee (1990) . Finally, the ability to predict damage in a composite 
laminate is demonstrated in section 4.4. 

4 . 1 Impact Response of a Steel Beam 

The impact response of a simply supported steel beam was predicted 
using both plane strain and brick elements to validate the transient 
solver. Beam, as shown in Figure 4.1, has dimensions of 10 x 1 x 1 in. 

A spherical steel impactor with equal mass to that of the beam, impacts 
the beam at its center with a velocity of 10 in/sec. The impactor is 
assumed to adhere to the beam during the impact. The following 
represents the material properties used for both impactor and the beam: 
Material : Steel, 

Density : 0.285 lbs/ in 3 , 

Poissons Ratio u : 


Youngs Modulus : 


0.29, 

30 x 10 6 lbs/in 2 . 
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Figure 4 . 1 

Simply Supported Steel Beam. 


As the model is symmetric, only half of the beam was modeled using 
both plane-strain and brick elements. A 2 x 2 Gaussian quadrature 
formula was used for plane strain elements and a 2 x 2 x 2 quadrature 
formula was used for the brick elements. The impact response of the 
model was obtained for an impact duration of 1000 |is . 

The deflection at the center of the beam is shown in Figure 4.2. 
The displacement values are about 5% lower than those obtained by 
Goldsmith (1960) , who replaced the beam with an equivalent spring-mass 
system. Since his formulation did not account for shear deformation in 
the beam, all impact energy was imparted to flexural mode greater 
deflections that should be observed in the actual impact. In the 
present case, both brick and plane strain elements include shear; hence 
some impact energy is channeled to the shear deformation mode, resulting 
in smaller deflections than that obtained by Goldsmith. 
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Figure 4 . 2 

Deflection at the Center of the Simply Supported Beam. 

Figure 4.2 also shows that the plane strain model yields smaller 
deflections than does the brick model. The reason for this discrepancy 
is that the plane strain assumption imposes an additional constraint 
that results in stiff er elements than the brick counterparts. 

A simple approximate analytical solution may be derived when the 
impactor mass is large compared to the mass of the beam. In such cases, 
the system is represented by a spring-mass system in which the mass of 
the beam constitutes little to the overall response and can be 
represented as one half the total beam mass lumped at the center. An 
equivalent spring with the stiffness equal to the flexural stiffness of 
the beam is used. The impactor mass is also lumped at the beam center. 


A half period for this model computed to be 955 jis whereas the finite 
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element model predicts 910 ns, with a difference of 5%. It should be 
recognized that spring-mass model yields an upper limit as no energy 
absorption due to shear deformation is assumed. On the other hand the 
finite element should form a lower limit due to the constraints imposed 
by the assumed shape functions . Higher-order elements should improve 
the accuracy of the finite element solution. Thus, the actual 
deflection and time period should not be further than 5% from the finite 


element solution. 
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4.2 T ransient Response of a Plate. 

To validate the composite element formulations, the transient 
responses of two types of plates was predicted. The first plate was 
made of an isotropic material . The second plate was made of a 
transversely isotropic laminated material. 

4.2.1 Xsotropic Plate 

A square plate, simply- supported, subjected to uniform pulse load 
is used to further validate the model . Due to symmetry, only a quarter- 
plate was modeled. The finite-element model, plate dimensions, and 
material properties used are shown in Figure 4.3, as follows: 


q(t) 

q 0 — 

SS : w=u=0 
CC: u=v=w=0 

t 

Dimensions =25X25X5 cm 
p = 8 X 10 * N sec /cm 4 
v = 0.25 

E = 2.1 X 10 N/cm 2 
, 10 N/SqCm 

4.3 

Finite Element Model of a Quarter Plate. 

(SS: Simply Supported; CC: Clamped) 
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A finite element model using brick elements with a mesh density of 
5x5x2 was used to model the quarter plate. The numerical 
computations were carried out using a 3 x 3 x 3 Gaussian quadrature 
formula. Figure 4.4, shows the deflection at the center of the plate. 
The deflection of the center obtained was very close to the result 
obtained by Reddy (1983) , even though the element formulations used were 
different. Figure 4.5, shows the variation of strain, kinetic energy, 
and the total system energy with time. 


I — — 

Displacement Vs Time 



Figure 4 . 4 

Transient Response at the Center of a Simply Supported Square 
Isotropic Plate subjected to Suddely Applied Uniform Pulse Load. 
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Figure 4 . 5 

Variation of Energies in a Square Isotropic Plate subjected to Suddenly 
Applied Uniform Pulse Load. (T.E : Total System Energy; K.E: Total 
System Kinetic Energy; STRN.E: Total System Strain Energy) 


When the load is applied, the kinetic energy increases as the 
nodal velocities increase. Simultaneously, the strain energy increases 
due to the resistance of the plate to the external force. The same 
resistance slows the nodal movements, causing a decrease in kinetic 
energy. Eventually, the strain energy becomes maximum when the kinetic 
energy vanishes. Afterwards, the nodes reverse direction and a mirror 


image is obtained. 
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4.2.2 Composite Plate 

The previous plate model was repeated using a composite material 
with E x = 25E 2 and G 12 = 0 . 5E 2 and E 2 = 2.1 x 10 6 N/cm 2 . Different lay-up 
sequences and two boundary conditions (clamped or simply supported) were 
tested: 

(i) f-45/451 — lay-up sequence: A finite element model with mesh 
density 5x5x2 was used to model a composite laminate with a [-45/45] 
lay-up sequence. Each ply along the thickness was represented by an 
element with the respective fiber orientation. A 3 x 3 x 3 Gaussian 
quadrature formula was used. Figure 4.6 as follows, shows the 
deflection at center of the plate for both boundary conditions : 



Figure 4 . 6 

Deflection at the Center of a Square Composite Plate with a Lay-Up 
Sequence [-45/45] subjected to Suddenly Applied Pressure Loading. 
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(ii) I3Q/45/90/P1 Lay-up Sequence; A mesh density of 3 x 3 x 1 

was used to model a laminate with [30/45/90/0] lay-up sequence. Also 
only one element was used along the laminate thickness direction. 

Figure 4.7, shows the center deflection of a clamped plate. 


Displacement Vs Time 



Figure 4 . 7 

Deflection at the Center of a Square Composite Plate with Lay-Up 
Sequence [30/45/90/0] subjected to Suddenly Applied Uniform Pulse Load. 


The results obtained from present software were close, and were 
comparable to those published by Reddy (1983) . The maximum difference 
observed was 6%. This variation was due to differences in element 
formulation, and the time integration schemes employed. These examples 
sufficiently demonstrate the accuracy of the present model. 
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4 - 3 I ne l ast i c I mpact R esponse of a Composite plate 

In the present work, inelastic impact was assumed because of its 
simplicity. However, Hertz contact law can be incorporated with some 
minor modifications to the software. A numerical case study from Lin 
and Lee (1990) was used to validate inelastic impact response of a 
composite plate. A clamped composite plate having dimensions 0.14 m x 
0.14 m with a lay-up sequence of [0 5 /90 5 /0 5 ] and the following material 
properties are used: 

E x = 40 Gpa, 

E 2 = 8.27 Gpa , 

^12 “ Gi3 “ G 23 = 4 ■ i 4 Gpa, 
v 12 = 0.26, 
h = 0.00335 Mt . , 
p = 1901.5 kg/m 3 . 

The impactor has a mass of 0.014175 kg and an initial velocity of 39.7 
m/s. Due to symmetry, only a quarter plate was modeled using a mesh 
density of 10 x 10 x 1. 2x2x1, and 2x2x2 integration schemes 

were used for computation of stiffness matrices of plies and element 
mass matrices, respectively. The displacement response, and the energy 
response, are shown in Figures 4.8 and 4.9, respectively. 
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Figure 4 . 8 

Deflection Response at the Center of a Clamped 
Composite Plate due to Inelastic Impact. 


Both the maximum deflection, and the time period, are lower than 
the results obtained by Lin and Lee (1990) by 10-15%. This difference 
can be attributed to the fact that the shell element formulation used by 
Lin and Lee does not include shear deformation. A better result may be 
obtained by using higher order approximation functions that culminate 
in a computationally larger model . 
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Figure 4 . 9 

Variation of Energies in a Clamped Composite Plate due to 
Inelastic Impact. (Strn.E: Strain Energy; 

K.E: Kinetic Energy; T.E: Total Energy) 

It is clear from Figure 4.9 that the total system energy is slowly 
increasing, though it is expected to remain constant (with no damping) . 
This increase in energy is due to numerical inaccuracies that are 
carried over from previous iterations. Accumulated numerical 
inaccuracies can cause the system to become unstable. The onset of 
instabilities can be observed in the strain energy and kinetic energy 
curves at the 142 |is mark. This unstable behavior can be avoided in two 
ways with some trade-offs: 

i) Use smaller time steps. However, a very small time step 
could result in computationally expensive analysis 

ii) Choose appropriate damping parameters to correct for the 
energy creepage as well as to account for structural 


damping . 
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4 - 4 Damage Prediction in a Composite Plate 

A model from Choi and Chang (1992) was analyzed in two case 
studies with the following conditions: 

i) Material degradation is not included in the analysis. 

ii) Material degradation is included in the analysis. Either 
complete or partial degradation could be applied. Complete 
degradation was assumed in this case study. 

4.4.1 No_Material De gradation 

A failure analysis of a model with no material degradation was 
solved. The problem consists of a spherical steel impactor moving with 
a velocity of 7.8 m/sec and stationary composite plate. The plate made 
of Fiberite T300/976 having the dimensions 10 cm x 7.6 cm x 0.403 cm and 
a lay-up sequence of [0 4 /-45 4 /45 4 /90 4 /45 4 /-45 4 /0 4 ] . The properties of the 
composite material are shown in Table 4.1. The impactor hits the plate 
at the center and is assumed to stick to the plate after impact. A 
finite element model of a quarter plate using a 5 x 4 x 1 size mesh was 
used. The model is shown in Figure 4.10. Stiffness computations were 
carried out using a 3 x 3 x 1 integration scheme while a 3 x 3 x 3 
integration scheme was used for computation of a diagonal mass matrix. A 
time step le-09 sec interval was used in this analysis. 
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Table 4 . 1 

Material Properties of Fiberite T300/976 Graphite /Epoxy 
source : (Choi and Chang 1992) . 
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Figure 4.10 

Quarter Plate Finite Element Model of a Clamped Composite Plate 
with Lay-Up Sequence of [0 4 /-45 4 /45 4 /90 4 /45 4 /-45 4 /0j . 


Displacement (mm) 
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Force Vs Time 



Figure 4 . 12 

Resisting Force Response of the Impact Point in Composite 
Plate Model with No Material Degradation 

The impact response obtained could not be verified as impact 
response for this particular example was not studied by Choi and Chang 
(1992) . However , the resisting force and deflection profiles obtained 
at the center of the plate due to impact were similar. 
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It is a well-known fact that damage in a system cannot occur 
without energy dissipation. The variation of energy with time is an 
important parameter to predict and to estimate damage. In the present 
case, the material damage was not included in the analysis, hence the 
total energy remained constant as expected as shown in Figure 4.13. 


Energy Vs Time 



Time ( V s) 


Figure 4.13 

Variation of Energies in a Composite Plate 
Model with No Material Degradation. 
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Layer 5 


[0 4 /45 4 /-45 4 /90 4 M5 4 /45 4 /0 4 ] 

1 


Layer 6 

[0/t5^5«/90«/-45 4 /45 4 /0 4 ] 

I 


(e) 


(f) 



(g) 

Figure 4.15 

Top View of the Predicted Damage with Model using 
No Material Degradation Method. 

Top View of Damage in (a) Layer 1 (b) Layer 2 (c) Layer 3 
(d) Layer 4 (e) Layer 5 (f) Layer 6 (g) All Layers 
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Time ( |i s) 


Figure 4 . 18 

Variation of Constituent Energies in a Composite 
Plate Model with Material Degradation. 


Energy dissipation caused by damage is illustrated in the 
Figure 4.18. The total energy curve indicates a gradual drop between 0 
and 100 |is, after which no further damage was incurred except at the 
very end. As there is no other source that consumes energy, damage is 
cause of the observed energy dissipation. Using Figure 4.18, the extent 
of damage may be estimated by the severity of the energy loss. Once the 
threshold value for energy loss is determined experimentally, the actual 
residual strength may be predicted from these curves. 
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4 . 5 Observations 

From the damage data obtained, and the sketches representing 
damage, the following observations were made: 

i) The severity of damage near vicinity of the impact is 
maximum. Almost all of the plies were damaged in this zone. 

ii) The damage seem to spread in the direction of fiber 
orientation, especially at the ply interfaces where the ply 
orientation changes. 

iii) Damage is most severe in bottom layers probably due to the 
higher bending stresses induced. It was also observed that 
the first point of damage is near the vicinity of the impact 
and in 0° and 45° layers at the bottom interface, i.e., 
between layer 1 and layer 2 . 

iv) The damage sketches indicate that there is no evidence of 
damage in the 7th layer of the no -degradation model, while 
the 7th layer of the degradation model exhibits the damage. 
The reason is that the material degradation model degrades 
the material when the damage is predicted thus softening the 
material . Due to the reduced stiffness in the damage 
layers, the unbalanced stresses are transferred to the 7th 
layer causing damage in this layer, as well. 

v) The damage predicted was more pervasive in the no material 
degradation model because no energy was dissipated during 
impact. In the material degradation model, material 
degradation resulted in energy dissipation, and caused the 


total system energy to be lowered. 
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CHAPTER 5 
Conclusions 


5 . 1 Conclusions 

An analytical model was developed to predict low-velocity impact 
damage in graphite/epoxy composite laminates. User friendly software 
package was developed based on the object-oriented implementation 
framework called Extensible Implementation Framework for Finite Elements 
(EIFFE) . 

Four case studies were analyzed first two case studies, i.e., 
simply supported steel beam and isotropic and transversely isotropic 
plate structures, validated the model. The third case study verified 
the inelastic response of the model. Damage prediction in 
graphite/epoxy composite laminates was studied in the fourth case study, 
which consists of two damage models. The first model did not include 
material degradation, while the second model included the material 
degradation. Overall, the predictions obtained from the model were 
comparable to results obtained by Choi and Chang (1992) . Based on the 
results obtained the following remarks can be made: 

i) Damage is severe in the vicinity of impact. 

ii) Damage is more prevalent in plies at the interface where the 
fiber orientation changes. 

iii) In-ply damage grows along the fiber orientation. From this, 
it can be construed that damage primarily occurs in the 
matrix and seem to grow in the fiber direction, causing 


delamination . 
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iv) The damage size predicted by the full material degradation 
model is less than the damage size predicted by the model 
with no material degradation. 


The model could be further improved to increase accuracy. 

Following are some of recommendations for improving the model : 

i) Incorporate Hertz contact model to more accurately model the 
impact forces. 

ii) Perform parametric study to determine the material 
degradation fraction and compare to experimental results. 

iii) Explore alternative material degradation methods; one such 
method is reducing material stiffness within the damaged 
elements during the analysis. 
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Program Listing 


// Module: COMPOSITE. H 

// Description: Interface file for Composite Element classes 

#ifndef COMPOS I TE_H 

#def ine COMPOS I TE_H 

#include <afx.h> 

#include <feobj .h> 

#include <struct rl . h> 

#include <material.h> 

#include <f stream. h> 

_CLASSDEF (RKC_3DCompositeTemplate) 

_CLASSDEF (RKC_Composi teElement ) 

_CLASSDEF (RKC_3DTemplate) 


class _CLASSTYPE RKC_3DTemplate : public VN_3DTemplate 

{ 

public : 

RKC_3DTemplate (const CString & Name, RVN_Material Material, 
PVN_FunctionArray pS FAr ray =NULL , 

PVN_Integrator pKIntegrator=NULL, PVN_Integrator 
pInitialFIntegrator=NULL, 

PVN_Integra tor pBodyFIntegra tor=NULL , PVN_Integrator 
pMIntegrator=NULL) ; 

virtual dataType Computes trainEnergy (RVN_StructuralElement 
rElement, RCVN_Vector rElementDisplacements) ; 

PVN_Vector GetCoordinatesAt (RVN_StructuralElement rElement, 
RCVN_DataArray X) ; 

protected : 

RKC_3DTemplate ( ) { } ; 

PVN_Vect or GetCoordinatesAt (RCVN_DataArray X) ; 


enum degrada t ionType { COMPLETE , PARTIAL , NOCOMPLETE } ; 

class _CLAS S TYPE RKC_3DCompositeTemplate : public RKC_3DTemplate 

{ 

DECLARE_ SERIAL (RKC_3DCompositeTemplate) 
public : 

RKC_3DCompositeTemplate (const CString & Name, RVN__Material 
Material , PVN_Integrator pLKIntegrator , PVN_Funct ionArray 
pS FAr ray =NULL , 
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PVN_Integrator pKIntegrator=NULL, PVN_Integrator 
pInitialFIntegrator=NULL, 

PVN_Integrator pBodyFIntegrator=NULL, PVN_Integrator 
pMIntegrator=NULL) ; 

~RKC_3DCompositeTemplate { ) ; 

RVNjSymMatrix GetE (dataType z) ; 

RVN_SymMatrix GetE ( ) ; 

indexType GetLayerNumAt (dataType z) const; 
indexType GetLayerNumAt (indexType PointNo) const ; 
virtual indexType NumLayerlntegPts () const ; 
virtual PVN_SymMatrix ComputeStif fnessMatrix (); 

PVN_SymMatrix ComputeStif fnessMatrixAt Point (indexType PointNo); 
virtual PVN_Vector ComputeAveragelnitialForceVector 
(RCVN_Vector ElementDisplacments) ; 
virtual PVNJVector ComputeAverageStressesAtPoint (RCVN_Vector 
ElementDisplacments , indexType LocalPointNo, indexType 
FrmLayerNo, indexType ToLayerNo) ; 
virtual PVNJVector ComputelnitialForceVector (RCVNJVector 
ElementDisplacements) ; 

PVN_Vector ComputeStressesAt (RCVNJVector ElementDisplacements # 
RCVN_DataArray X) ; 

void GetlntrinsicCoordinatesAt (RVNJDataArray X, indexType 
PointNo) const; 

BOOL UpdateElementCondition (RCVN_Vector 

ElementDisplacments , degradationType eDegrade) ; 
void Serialize (CArchive &ar); 

protected: 

dataType prevOrientation, m_dStrainEnergy; 

PVN_SymMatrix m_pEp; 

PVN_NdxArray m_pMarks Array; 

PVN_Matrix m_pGaussPtsWtsMatrix; 

PVN_Cube Integrator mjpLKIntegrator ; 
void CreateGaussPtsWtsMatrix ( ) ; 
void CreateMarksArray ( ) ; 
virtual void ComputeEO; 

RKC_3DCompositeTemplate ( ) ; 

RKC_3DCompositeTemplate (RCRKC_3DCompositeTemplate ) ; 
RCRKC__3DCompositeTemplate operator= (RCRKC_3DCompositeTemplate 

) ; 

}; 


class _CLASSTYPE RKC_CompositeElement : public VN_StructuralElement 

{ 

DECLARE__SERIAL (RKC_CompositeElement ) 
public : 

RKC_CompositeElement (RRKC_3DCompositeTemplate 
Template , PVN_NodeArr ay pNodeArray, 

RVN_DataArray rThicknessArray , RVN_DataArray 
rOrientationArray) ; 

~RKC_CompositeElement () ; 
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indexType NumLayersO const; 
dataType Thickness ( ) const ; 

dataType LayerThicknessAt (indexType LayerNo) const; 
dataType LayerOrientationAt (indexType LayerNo) const; 
dataType LayerLowerCoordinateAt (indexType LayerNo) const ; 
dataType LayerUpperCoordinateAt (indexType LayerNo) const ; 
dataType LayerThicknessRatioAt (indexType LayerNo) const; 
dataType LayerMidZCoordinateAt (indexType LayerNo) const ; 
dataType LayerMidEtaCoordinateAt (indexType LayerNo) const; 
indexType GetLayerNumAt (dataType z) const; 
dataType GetlntrinsicCoordinate (dataType z) const; 
dataType GetExtrinsicCoordinate (dataType Eta) const; 
dataType Transf ormCoordinate (dataType CoordX, dataType 

LowerLimitX, dataType UpperLimitx, dataType LowerLimitx, 
dataType UpperLimitx) const; 

PVN_SymMatrix Computes t if fnessMatrixAtPoint (indexType PointNo) 
const ; 

void ComputeAveragelnitialForceVector (RCVN_Vector 
ElementVector) ; 

RVN_Vector GetAveragelnitialForceVector ( ) const ; 
virtual void ComputelnitialForceVector (RCVN_Vector 
ElementVector) ; 

void RemoveAveragelnitialForceVector ( ) ; 

BOOL UpdateElementCondition (RCVN_Vector ElementDisplacements , 
degradationType eDegrade = NOCOMPLETE) ; 
void UpdateStif fnessMatrix (dataType DegradationFactor = 1.0); 
BOOL GaussPtStsOK ( indexType PointNo) const; 
void SwitchAllGaussPtsSts (BOOL bFlag); 
void SwitchGaussPtSts (indexType PointNo, BOOL bFlag); 
indexType GaussPtReport (indexType PointNo) const; 
void SwitchGaussPtRptSts (indexType PointNo, indexType Flag) ; 
void SwitchAllGaussPtsRptSts (indexType Flag); 

RVN_SymMatrix GetE ( ) const ; 

virtual RVN_SymMatrix GetE (dataType z) const; 
void PrintElementCondition (of stream &os, indexType 
ElementNo, indexType iter) ; 
void Serialize (CArchive &ar) ; 
protected : 

indexType m_uNumLayers , m_iKIntegPts , m_iLKIntegPt s ; 
dataType m_dThickness , m_dLowerLimit , m_dUpperLimi t ; 
PVN_DataArray m_pThickArray, m_pOrientationArray ; 

PVN_NdxArray m_pbGaussPtsPrntStsArray ; 

BOOL *m_pbGaussPtsStsArray ; 
dataType ComputeThickness ( ) ; 
void CreateGaussPtsStsArray ( ) ; 

RKC_CompositeElement { ) ; 

RKC_CompositeElement (RCRKC_CompositeElement ) ; 
RCRKC_CompositeElement operator= (RCRKC_CompositeElement ) ; 

} ; 

#include "compos it . ini" 

#endif 
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// 

// 

// 

// 

n 


Module: INITIALC.H 

Description: Interface file for Intial 

Boundary Condition class 


#ifndef _INITIALC H 

#define _INITIALC H 

#include <afx.h> 

#include <datadefs.h> 

#include <matrix.h> 

^include <feobj.h> 

_CLASSDEF (RKC_InitialC) 

class _CLAS STYPE RKC_InitialC : public CObject 

{ 

public : 

RKC_InitialC (RVN_Node rNode, dataType Value, UINT dir); 
RVN_Node NodeO const; 
dataType Value () const; 

UINT Direction {) const; 

PRKC_InitialC CopyTo (RVN_Node newNode) ; 
void Apply (RVN_Vector F) const; 
virtual void pr intOn (ostream Sc os) const; 
protected : 

RKC_InitialC ( ) ; 

PVN_Node m_j)Node ; 
dataType m_dValue; 

UINT m_iDirection; 
private : 

RKC_InitialC (RCRKC_InitialC) ; 

RCRKC_InitialC operator = (RCRKC_InitialC) ; 

}; 

#include " InitialC . ini M 
#endif 



Chenna 83 


// Module: 

// Description 

// 

#ifndef TRNSCOLL_H 

#define TRNSCOLL_H 

#include <f ecolls . h> 

#include "initialc . h" 

_CLASSDEF (RKC_InitialCArray) 

_CLASSDEF (RKC_InitialCList ) 

class _CLASSTYPE RKC_InitialCArray : public VN_ObjectArray 
public : 

RKC_InitialCArray (BOOL bOwnElements = FALSE) : VN_Object Array 
(bOwnElement s ) { } ; 

RRKC_InitialC IC(int i) { return (RRKC_InitialC) * Element At (i) ; 

} 

int AddIC (PRKC_InitialC pIC) { return Add ( (CObject *) pic) ; } 

void SetICAt (int i, PRKC_InitialC pIC) { SetAt (i, (CObject *) 
PIC) ; } 

void InsertICAt (int i, PRKC_InitialC pIC, int nCount = l) 

{ InsertAt (i, (CObject *) pIC, nCount); } 

void SetICAtGrow (int i, PRKC_InitialC pIC) { SetAtGrow (i, 
(CObject *) pIC) ; } 

}; 

class _CLASSTYPE RKC_InitialCList : public VN ObjectList 

{ 

public : 

RKC_InitialCList (BOOL bOwnElements = TRUE) : VN_ObjectList 
(bOwnElements) { } ; 

POSITION AddHeadIC (PRKC_InitialC pIC) { return AddHead 
((CObject *) pIC) ; } 

POSITION AddTaillC (PRKC_InitialC pIC) { return AddTail 
{ (CObject *) pIC) ; } 

POSITION InsertICAfter (POSITION pos, PRKC_InitialC pIC) { 
return InsertAfter (pos, (CObject *) pIC) ; } 

POSITION InsertICBefore (POSITION pos, PRKC_InitialC pIC) { 
return InsertBefore (pos, (CObject *) pIC) ,- } 
void SetICAt (POSITION pos, PRKC_InitialC pIC) { SetAt (pos, 
(CObject *) pIC) ; } 

PRKC_InitialC RemoveHeadIC ( ) { return (PRKC_InitialC) 

RemoveHead ( ) ; } 

PRKC_InitialC RemoveTaillC ( ) { return (PRKC_InitialC) 

RemoveTail ( ) ; } 

RRKC_InitialC IC (POSITION pos) { return (RRKC_InitialC) 

♦GetAt (pos) ; } 

RRKC_InitialC HeadICO { return (RRKC_InitialC) *GetHead(); } 
RRKC_InitialC TaillCO { return (RRKC_InitialC) *GetTail(); } 


TRNSCOLLS . H 

Interface file for Transient 
Boundary Condition Collection classes 
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RRKC_In.it ialC NextIC (POSITION & pos) { return (RRKC_InitialC) 
*GetNext (pos) ; } 

RRKC_InitialC PrevIC (POSITION & pos) { return (RRKC_InitialC) 
*GetPrev (pos) ; } 

}; 

#endif 


// Module: 

// Description 

#ifndef TRANS SOL_H 

#def ine TRANS SOL_H 

#include <afx.h> 

#include <feobj .h> 

#include <structrl . h> 

#include <material . h> 

#include" trnscoll .h" 

#include "composit . h" 

#include <f stream. h> 

_CLASSDEF (RKC_TransientModel) 

enum approachType { DIRECTVECTOR , DEFAULT } ; 
enum force Type {AVERAGE FORCES , DEFAULTFORCES } ; 

class _CLASSTYPE RKC_TransientModel : public VN_StructuralModel 

{ 

public : 

RKC_TransientModel (const CString & Name=" " ) ; 

~RKC_TransientModel () ; 

virtual void FlushO; 

void AddlnitialC (PRKC_InitialC pIC) ; 

int NumlnitialCs () const ; 

void ComputeAverageElementlnitialForces () ; 
virtual void Ini tiate Impulse (da taType dMass , approachType 
eAppType = DEFAULT, for ceType eForceType = 

DEFAULTFORCES , dataType Alpha =0, dataType Beta = 0) ; 
virtual void Solve (dataType Incr) ; 
virtual void ComputeNodalAccelarations (RVN_Vector 

rAccelarations , RCVN_Vector InvM, RCVN_Vector F, RCVN_Vector 
IVel ) ; 

virtual void ComputeNodalAccelarations (RVN_Vector 
rAccelarations, RCVN_Vector InvM, RCVN_Vector F) ; 
virtual void UpdateNodalDisplacements (RVN_Vector 

rSolution, RCVN_Vector IVel, RCVN_Vector Acc, dataType 
DeltaTm) ; 

virtual void UpdateNodalVelocities (RVN_Vector rVelocites, 
RCVN_Vector rAcc, dataType DeltaTm) ; 


TRANSSOL. H 

Interface file for Transient Model class 
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dataType ComputeCriticalTime (dataType TOL = le-7, indexType 
Iter = 2000) ; 

dataType ComputeSystemKineticEnergy () ; 
dataType ComputeSystemStrainEnergy () ; 

RVN_Vector Nodal Velocities { ) const ; 

RVN_Vector Nodal Forces ( ) const ; 

RVN__Vector NodalAccelarations ( ) const ; 
virtual void UpdateModel ( ) ; 

BOOL UpdateElementsCondition (degradationType eDegrade = 
NOCOMPLETE , dataType DegradationFactor=l) ; 
void SerializeVectors (CArchive & ar) ; 

protected: 

RKC_InitialCList m_InitialCList ; 

PVN_Vector m _pIVel , m_pAcc, m_pInvMass, 

m _pGlobal Ini tialForceVector , m_pGlobalKXVector , 
m_pGlobalCXVector , m_pGlobalForceVector , m_pMassVector ; 
approachType m_eAppType ; 
forceType m_eForceType ; 
dataType m_dAlpha / m_dBeta; 
of stream os; 

BOOL m_bCrTmComputed; 

void ApplyEssentialBCs (RVN_Vector x, applyType Flag) ; 
void AssembleMassVector ( ) ; 

void ApplylnitialConditions (RVN_Vector Value, dataType Mass); 
void ApplyNaturalBCs (RVN__Vector Value); 
void AssembleElementlnitialForces () ; 
void ComputeAssembleElementKXVectors (RCVN_Vector 
ElementDisplacments) ; 

void ComputeAssembleElementCXVectors (RCVN_Vector 
ElementVelocityVector ) ; 
void UpdateElementlnitialForces () ; 

void UpdateElementStif fnesses (dataType DegradationFactor=l) ; 

void UpdateForceVector () ; 

void InvertMassVector ( ) ; 

void CreateVectors ( ) ; 

void RemoveVectors ( ) ; 

void RemoveElementStif fnesses ( ) ; 

private : 

RKC_Trans ientModel (RCRKC_TransientModel) ; 

RCRKC_Trans ientModel operator= (RCRKC_TransientModel ) ; 


#include "transsol.ini 
#endif 
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// Module: MODEL. H 

// Description: Interface file for Conversion Module From 

// Nastran Data File to EIFFE Objects 

#ifndef _MODEL_H 
# define MODEL H 


#include 

#include 

#include 

^include 

#include 

#include 

#include 

#include 

#include 

#include 

#include 

#include 

#include 

#include 


<afx.h> 
<afxwin . h> 
<ostream. h> 

< f stream. h> 
<stdio . h> 
<stdlib. h> 
<string . h> 
<iomanip . h> 
<f eobj . h> 
<structrl .h> 
"nasdef .h M 
"compos it . h" 
"transsol . h" 
"plancomp . h" 


_CLASSDEF (RKC_Model) 

class _CLASSTYPE RKC_Model 

{ 

public : 

RKC_Model (model Type eModelType, solutionType eSolType, CString 
ProblemHeader) ; 

~RKC_Model ( ) ; 

PVN_ElementTemplate CreateElementTemplate (modelType eModelType, 
elementType eElementType , materialType eMat ) ; 
PVN_StructuralTemplate 

Creates tructuralElementTemplate (elementType eElementType, 
materialType eMat) ; 

PVN_Element CreateElement (RVN_Element Template pElementTemplate , 
PVN_NodeArray pConNode Array, materialType eMat) ; 
PVN_StructuralElement 

CreateStructuralElement ( RVN_S t rue tur a 1 Template 
pElementTemplate , PVN_NodeArray pConNodeArray , materialType 
eMat) ; 

void SolveStaticModel ( ) ; 
void SolveTransientModel ( ) ; 

void ReadinNastranModel (if stream _FAR & os, CString & 

Mat Filename) ; 

void PrintNodal (RVN_Vector rVector, CString Item,ostream _FAR 
&OS , indexType Case = 1 , print FormatType 
eFormat=SPECIAL) const ; 

void PrintStresses (ostream _FAR &os, indexType Case= 1, 
printFormatType eFormatType = SPECITUj) const; 

RVN_Vector Solution 0; 
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void GetViewClass (CView *pView) ; 
void ComputeCriticalTimeStep ( ) ; 

dataType ComputeElementStrainEnergy (indexType ElementNo) ; 

protected : 

CString m_sName; 
indexType m_uNumI t era t ions ; 

UINT m_iDOFNum; 

PVN_Vector m__pSolut ionVector , m_pVelVector , m_jpAccVector , 
m_pForVector ; 

CView *m_pView; 

PVN_FunctionArray m_pFunctionArray ; 

PVN_Integrator m_pKIntegrator , m_pInitialFIntegrator , 
m _j>BodyFIntegra tor , m_pMass Integrator , m_pLKInt egra tor ; 
PVN_StructuralModel m__pFEModel; 

PVN_Element Temp late m_pElement Template ; 

PVN_Material m_pMaterial ; 
indexType m_iPrevElement ; 

indexType m_iElementID, m_iNodeNo, m_iStressElementID, m_iLayerNo; 
dataType m_dDel taTrn, 

m_dCrTm, m_dStrainEnergy , m_dKineticEnergy , m_dModelTime ; 
modelType m_eModelType ; 
solutionType m_eSol ut ion Type ; 
elementType m_e Element Type ; 
materialType m_eMaterialType ; 

PVN_NodeArray m_pNodeArray; 

PVN__Element Array m_jpEl erne nt Array; 

PVN_ElementTemplateArray m_pElementTemplateArray ; 

PVN_Ob j ect Array 

m_pNdxContainer , m_pOrDataContainer # m_j?ThDataContainer ; 
PVN_NdxArray m jpS t res sElement Array; 
of stream 

osdisp, osvel , os force, osacc, os ini # osdam, oscomb, osstress ; 
void FirstScan (if stream _FAR & is); 
virtual void ReadGrid (if stream _FAR & is); 
virtual void ReadElement ( if stream _FAR & is, CString & 

Mat Filename) ; 

virtual void SetSizeof Arrays () ; 

virtual void ReadBoun da ryCondit ions (if stream _FAR 
&is , boundaryType eBoundary) ; 

virtual void ReadElementProperties (const CString & sMaterial , 
CString & MatFilename) ; 

virtual PVN_Element Creates he 11 El ement (RVN_UINTAr ray iConnData, 
materialType eMat) ; 

virtual PVN_Element Create8NodeBrickElement (RVN_UINTArray 
iConnData, materialType eMat); 

virtual PVN_Integrator CreateKIntegrator (elementType eElement); 
PVN_Integrator CreatelnitialBodyFIntegrator (elementType 
eElement) ; 

PVN_Integrator Great eBodyFIntegrator (elementType eElement); 
PVN_Integrator CreateMass Integrator (elementType eElement); 
PVN_Integrator CreateLayerKIntegrator (elementType eElement) ; 
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PVN_FunctionArray CreateShapeFunctionArray (elementType 
eElement) ; 

void Createlntegrators (elementType eElement); 
virtual PVN_3DTemplate Create3DTemplate (mater ialType eMat); 
virtual PVN_PlaneStressTemplate CreatePlaneStressTemplate 
(mater ialType eMat) ; 

const char* Processstring (CSt ring _FAR &var) ; 

virtual void ReadLabel ( is t ream _FAR & is, const char* Label); 
void CreateOutputFiles ( ) ; 
void CloseOutputFiles ( ) ; 

void PrintAllVectors (indexType Everylter=l) ; 

void Print (indexType NodeNo, RVN_Vector Displacements, 

RVN_Vector Velocities, RVN_Vector Accelarations , RVN_Vector 
Forces , of stream _FAR & os, indexType Case ) const; 
void PrintStressesAt (indexType ElementNo, RCVN_Vector X, 

RCVNJVector rStressVector , RCVN_Vector rStrainVector , of stream 
_FAR & os, indexType Case) const ; 
void PrintElementsStatus (of stream &os, indexType iter); 
PVN__Vector ComputeElementStressesIn (RVN_StructuralElement 
rElement, RVN_DataArray X, indexType LayerNo=0) ; 

} ; 

#endif 


c-sl 




